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Abstract 

Background: We aimed to investigate the effect of cell-cell dipole interactions in the 
equilibrium distributions in dielectrophoretic devices. 

Methods: We used a three dimensional coupled Monte Carlo-Poisson method to 
theoretically study the final distribution of a system of uncharged polarizable particles 
suspended in a static liquid medium under the action of an oscillating non-uniform 
electric field generated by polynomial electrodes. The simulated distributions have 
been compared with experimental ones observed in the case of MDA-MB-231 cells 
in the same operating conditions. 

Results: The real and simulated distributions are consistent. In both cases the cells 
distribution near the electrodes is dominated by cell-cell dipole interactions which 
generate long chains. 

Conclusions: The agreement between real and simulated cells' distributions 
demonstrate the method's reliability. The distribution are dominated by cell-cell 
dipole interactions even at low density regimes (10 5 cell/ml). An improved estimate 
for the density threshold governing the interaction free regime is suggested. 
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Background 

As first defined by Pohl [1,2], the term " dielectrophoresis" is used to describe the 
"ponderomotive" force exerted by a non-uniform electric field on polarizable neutral 
particles. Such force allows for the controlled manipulation of micro and nano-sized 
particles dispersed in colloidal solutions. Application fields include: cell partitioning 
and isolation [3,4], bio-structure assembling [5], nanostructure (e.g. carbon nanotube) 
deposition [6], filtration systems for oils purification [7] etc. Among these, the separation 
of rare cells [8] is a specifically promising one as dielectrophoresis allows the capture/ 
separations of cells without the use of biomarkers; relying, instead, in the strong 
selectivity of the dielectrophoretic (DEP) response [9] which depends on the particle 
mass, shape and composition. Indeed, recently, this selectivity has permitted to discrimin- 
ate the tumor cell types of the NCI-60 panel from Peripheral Blood MonoNuclear cells 
(PBMNs) [10]. However, although many intriguing micro-structures have been fabricated 
in research laboratories, DEP devices have hardly gone beyond the proof-of-concept stage 
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[11]. One of the problems that are hindering development and engineerization of the 
devices is the limited use of accurate numerical tools for their design which, in turn, is 
due to the computational complications arising by the particle-particle dipole interaction. 
Indeed, particle kinetics (i.e. the particle velocity field) in DEP devices can be easily 
calculated by mean of Poisson solvers and direct integration of the equation of motion 
only in the non-interacting particle approximation [12]. This approximation is not 
valid in the accumulation regions of the DEP devices where, due to the increased 
particle concentration, dipole-dipole interactions become important and can promote 
the formation of clusters and significant rearrangements of the particle space distribution 
[13-16]. These many-particle effects can be accurately simulated solving directly the 
equations of motion in the few-particles limit [17] i.e. this approach is not applicable 
for the simulation and design of realistic systems. Another possible approach is the use 
of reaction-diffusion models [18-20] but this approach needs an "ad hoc" parameter 
calibration to effectively consider the dipole-dipole interactions in compact models. 

Recently a coupled Monte Carlo-Poisson (MC-P) method [21] has been implemented 
which allows simulating a large number of particles in large active zones (within the 
experimental range), explicitly including particle-particle interactions. The MC-P 
method has pointed out the relevance of this inclusion in the modeling predictions for 



the simplified condition of Two Dimensional (2D) electric field E [~r\ distribution, 



where E (j\ explicitly depends on two space coordinates E ( ~r ] = E (x,y) as in the case 



of very long interdigitated electrodes. However, the possibility to apply this approach 
for the numerical design of devices exploiting more complex fully Three Dimensional 
(3D) electric field distributions has not been yet demonstrated. Moreover the MC-P 
predictions have never been compared with real cell distributions in dielectrophoretic 
devices, in order to confirm their reliability. Aiming to the two objectives of the model 
extension and validation, we have improved the application potentiality of the MC-P 
method to simulate the features of devices generating 3D electric field distributions. In 
addition we applied the simulation method to the case of polynomial electrodes which 
are known to produce well defined 3D non-uniform electric fields and are used for the 
study of negative dielectrophoresis [22] or for the determination of particle dielectrophoretic 
response through electrorotation analysis [23,24]. We compare the simulated results with 
experimental distributions obtained in the same electrodes geometry to evaluate the role of 
p-p interactions and definitively demonstrate the predictive potential of this methodology. 

Methods 

Computational algorithm 

A detailed description of the method can be found in ref. [21], here we summarize the 
key aspects of the simulations, specifically focusing on the 3D implementation. 
In the first-order approximation of the polarization, an isolated particle immersed in a 

media and subjected to a non-uniform electric field = E (j^j exp(/a)£) oscillating at 

the /= co/2tt frequency, will be subjected to an effective averaged potential energy [25]: 






(i) 
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where E rms (j^j = £ max {l^j / y/2 is the effective value of the varying electric field and 
a e ff is the average polarizability of the particle defined as: 



a e ff = 3VRe(s m )Re 



El-E„ 



si + 2e n 



3VRe(e m )Re\f CM (oj)} (2) 



where Vis the particle volume, = s^ m -ia^ m /co are the complex dielectric constants of 
the particle (1) and the media (m) and fcM(°>) is the Clausius-Mossotti factor which fully 
characterizes the dielectric response of the particle in the given medium. The isolated par- 
ticle approximation holds only in the diluted density limit, i.e. only if the average distance 
between two particles in the colloidal solution is always very large otherwise an effective 
particle-particle interaction has to be considered as a result of the local distortion of the 
electric field lines generated by itself and by the other particles. In this case, the DEP force 
acting on each suspended particle can be directly calculated by means of the Maxwell ten- 
sor [9]: T=Tij = s m (EfEf*- 0.5<%£f £f ) over the closing surface of the particles: 



t-DEP 



■ j^T-n^jdA (3) 



Note that whereas the E field in Eq. 1 is the field generated by the external electrodes 

only, E in the Eq. 3 must be calculated considering all the particle presence. This dir- 
ect calculation is not practically feasible in the kinetic simulation of large systems, since 
the particle distribution continuously changes in the space requiring an integration 

of the Poisson equation, WE tot = 0, at each simulation step. A more efficient ap- 
proach, that requires only the evaluation of the external field, can be implemented 
approximating the total distorted electric field with the sum of the field generated by 
the external electrodes plus the contributions of the dipoles induced in all the particles: 

all the particles 

e"(t5)-e(t;)+ £ ^(r,) [26]. 

The reliability of this approximation has been demonstrated in Ref. [13] with the aid 
of the full calculation based on Eq. 3 for the case of two spherical particles immersed 
in a uniform external electric field: the magnitude, the angular dependency and the 
scaling with the distance of the calculated force are similar to those derived in the 
interacting dipoles approximation [26]. 

The effective potential energy that generalizes Eq. 1 for the case of particle-particle 
instantaneous interactions, in the dipole approximation when the multipole terms and 
the mutual polarization of the particles can be neglected, is: 

U ir ~-\Ke( Vi (7).^)*) =-iRe(p / (7)).E i (7))*) (4) 

where E ; (j^j (e, (j^j ) is the electrical field generated by the dipole of the particle j (i) 
at the position > while P/ (~^) (p/ (o) ) * s tne ^ip°l e moment induced on the 



i-eme (j-eme) particle by the external field E^j(E(r ; )) . The electric field in r- 
generated by a particle in r y is: 



3 n [n • p.j-p. 

3 



47rRe(£ m ) R t) 



(5) 
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where n— Ry/Rij which leads to the following expression for the average effective 
potential energy: 

1 1-3 cos ( 0\ j cosf ff . J . v _^ f \\ 

where ^ , are the angles between the vectors E (j^ji E {j^j > an d n an d ^? a ig 
are the average polarizations for the i and / particles. 

We used a Monte Carlo methodology to efficiently determine the equilibration of a 
3D system of interacting particles suspended in a static liquid medium under the action 
of an oscillating non-uniform electric field generated by polynomial electrodes. The 
particles are considered as hard spheres with radius r t and the configuration energy is 

£({n, ...r„}) = J2 U eff (3) +E Q « ft' 7 !) ( ? ) 

i ij 

where U e ff and Uy are calculated by means of the equations (1) and (6). 

The solution of the Poisson equation, which allows determining the 3D electric field 
spatial distribution, has been calculated by means of a finite element numerical solver 
[27]. Then, the numerical estimate of the electric field has been interpolated in the 
cubic grid of the MC simulation box and, in order to obtain an accurate resolution, the 
distance of next-neighbor nodes in the simulation grid has been set equal to the particle 
radius. Equilibration kinetics, which is reached by a sequence of stochastic jumps for the 
simulated particles, can be inferred from the simulation if Brownian motion and 
particle-particle scattering due to the hard-sphere behavior lead to an approximated 
diffusive behavior. In this case we can estimate the time interval At between two 
consecutive displacement events as 

NAt = 1/6 x (Ad) 2 /D - 1/6 x (Ad) 2 /D Brow = 1/6 x (Ad) 2 (3nrja)/k B T (8) 

where D « D Bw w is the effective diffusivity, N is the number of simulated particles, Ad is 
the elementary particle displacement, rj is the medium viscosity, a is the particle radius, 
k B is the Boltzmann constant and T is the system temperature. 

Experimental setup 

The polynomial electrode design described in the previous section, has been fabricated 
by deposition of 10 nm of Titanium followed by 200 nm of Nickel on a standard micro- 
scope glass. The electrodes were delineated by lithographic methods followed by wet 
etching. The device has been energized using a Protek 9205C signal generator which 
applied, consistently with the simulated systems, a sinusoidal voltage signal of 8V pp 
value at 1 MHz for 180 sec (long time allow for cells equilibration). The final distribution 
was observed with a standard lOx phase contrast inverted microscope. The human breast 
cancer cell line MDA-MB-231 were cultured according to American Type Culture 
Collection (ATCC) instructions. The cells, just before DEP tests, were suspended in 
a low conductive buffer (used as the elute in all our experiments) composed of 9.5% 
ultrapure sucrose (S7903, Sigma- Aldrich), 0.3% dextrose (Fisher D-16), and 0.1% 
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Pluronic F68 (P1300, Sigma-Aldrich) titrated to a conductivity of 30 mS/m (consistent 
with Monte Carlo simulations) by adding KG with the aid of a conductivity meter. The 
buffer had an osmolarity of 320 mOs/L and a pH of 7 and the experiments were con- 
ducted at room temperature (22°C). The cells, suspended in DEP buffer at concentration 
of 5xl0 5 cells/ml were pipetted into the chamber and occupied a total volume of about 
100 \A when a cover slip was placed over the rubber o-ring. 



Results 

In order to study the interplay between the dielectrophoresis response and the particle- 
particle dipole interaction in a realistic case, we consider a suitable dielectric model of a 
colloidal system of MDA-MB-231 cells dispersed in a weakly conducting water solution 
(see later). The complex dielectric constant s e ff of these cells can be approximated by the 
solution of the following system of equations: 

_ _ 1 + 2dr 3 A 
£eff ~ 1-drU 

(9) 

dr = A = — 

where £/« jme = £in,me-i G in,mel°) are the inner cell and membrane complex dielectric 
constants and a = 6.2 \xm and d- 10 nm are the cell radius and membrane thickness. 
The following values have been used: s me = 24£ 0 , £ in = 50£o, &m = 30 rnS/rn, Oi n = 
0.2 S/m 1 a me = 10" 7 S/rn, [28], and the calculated real part of the factor /cm(^) is shown 
in Figure 1. 




10 



10 



10 10 10 

frequency (hz) 



10' 



10° 



Figure 1 Real part of the Clausius-Mosotti factor /bwM calculated with the dielectric model of the 
MDA-MB-231 cell at a medium conductivity of 30 mS/m. At the frequency of f= oo/2n = 10 6 Hz, used 
for both simulations and experiments, the cells are subjected to a positive-DEP, i.e. they will tend to move 
towards high electric field regions. 
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The cells are free to move in a cubic computational box of dimensions 
1600x1600x1500 um 3 with four polynomial electrodes located at the bottom of the 
box (see Figure 2) whose shapes can be described by the following parametric system: 



D<x<L 



±Vx 2 +D 2 



(10) 



for each electrode [22]. D represents half the distance of opposing electrodes whereas 
L is related to the electrode width. Referring to Figure 2, right and using Eq. 10, we 

have that L\ = V2D and L 2 = v^-V^-D 2 ). To improve the capturing efficiency of 

the DEP system the electric fields in p 2 and p 2 must be of the same order (otherwise 
the trapping regions will be limited to the p 2 regions). Given that |£rms|(Pi) * V PP IL\ 
and \E rms \ (p 2 ) * Vpp/Li D must be of the same order of L, in this study we sat D = 390 urn 
and L = 460 um (note that in ref. [22] D = 64 um and L = 90 um which favored negative- 
DEP only). Figure 2, left shows the simulated box with the four electrodes situated at the 

bottom. The intensity map of \E rms \ (j^j is also reported in grey scale for the oscillating 

four-electrodes configuration at a frequency of 1 MHz and a peak value of 8 V for V pp 
with 180° phase difference between neighbor electrodes (i.e. V\(t) = - V 2 (t) = V pp sin(2nft)). 
Figure 2, right shows the electric field at 50 um from the bottom surface, associated to the 
considered system. The regions of high electric field (to which the cells will move under 
positive-DEP) are located all around the edges of the electrodes. In the considered system 
the highest electrical fields are located at^ 7 and p 2 equivalent points. 

In Figure 3 we show the initial random distribution of 1920 cells (corresponding to 
a density of 5xl0 5 cells/ml) (left) and the final (right) equilibrium condition after 
2xl0 8 Monte Carlo iterations. The final distribution is the result of the minimization 
of Eq. 7, which induce a movement towards the p 2 and p 2 regions (minimization of 
the first term) and an aligned of the cells (causing cell-chains) along the electric field 
direction (minimization of the second term). Specifically particles tend to form 




Figure 2 Left, side view of the simulated computational box (1600x1600x1500 um 3 ) with, 
superimpose in grayscale, intensity map of |E rms | (7^j. Right, Top view of the system, the intensity map 

of |E rms | (7^J shows the areas of higher electric field (brighter areas). L 1 and L 2 are the distances between 
the opposing electrodes in p 1 and p 2l respectively. 
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. . t 



Figure 3 Initial (/eft) and final (r/g/if) spatial distribution of N = 1920 MDA-MB-231 cells 
suspended in a saline water solution of conductivity of 30mS/m subjected to an oscillating 
electric field f"=10 6 Hz. The cell distribution is obtained starting from a random distribution after 
2x70 s MC iterations. Dipole chains are evident at the bottom of the simulation box and are 
associated to high regions of electric field. 



larger chains in the regions of more intense electric field which are located mainly 
around p 2 equivalent points, similar distributions have been observed for densities as 
low as lxlO 5 cells/ml (not shown). It is important to note that this tendency does 
not depends on the DEP polarity because the a l e ^d e ^ term in Eq. 7 is always positive 
for identical cells (whereas the a l e g term in Eq. 1 can be positive or negative, depending 

on the value of fcM^))* 

From these results we can infer that particle-particle interactions compete with 
the dielectrophoretic force-field, which would otherwise massively trap (in p-DEP 
conditions) the particles in the regions where the gradient of the electric field is 
larger. Note also the cells in the region far away from the electrodes which are not 
trapped by the DEP field, this allows for the definition of a depletion volume as the 
region where cells are effectively attracted to the electrodes. The connection between 
depletion region, electrodes geometry and particle-particle interaction is currently 
under investigation. 

Figure 4 shows a comparison between the simulation and the experimental cells 
distribution. As can be seen the two distributions are equivalent when the statistical 
approach of the equilibration is considered. In both cases the cells distribution near the 
electrodes is dominated by cell-cell dipole interactions which generate long chains. 
Note also the, relative small, discrepancy between the simulated and experimental cell 
distribution in the region at the center of the four electrodes. One possible explanation 
of this weak discrepancy is the current preliminary calibration of the friction forces 
between the cells and the wall on the bottom of the chamber (this force is included in 
the MC code). Work is underway to improve the parameter calibration and definitively 
refine the simulation results. The highest value of the cell concentration to avoid 
particle-particle interactions strongly depends on the electrodes and system geometry, 
on the particles dimensions and on the polarization factors. As a consequence a general 
prescription to neglect dipole-dipole interaction in the design of a device cannot be 
easily given. Jones [29] showed that cell- cell interactions should become significant 
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Figure 4 left) Top view of the final spatial distribution of the Monte Carlo simulation (same of 

Figure 2) with superimposed the intensity map of |E rms | (7^j (same as Figure 2, right). The association 

between dipole chains and high regions of electric field is evident. Right) Top view of experimental spatial 
distribution of MDA-MB-231 after 180 sec of DEP manipulation. Black areas are associated to the metal 
electrodes. The chains are clearly present all around the electrodes. 



when the ratio of cell spacing 5 to cell diameter d falls below ~5. In a simple approxi- 
mation we could assume that all cells will eventually settle to the bottom of the 
chamber, which gives: 



S/d = l/(dh l ' 2 C 1 ^ >5^C< C threashold 



- j^r cells/ml 



(ii) 



where h is the chamber height and C is the cell concentration. Setting h = 1500 \im 
and d = 2xr= YlA\Avn we obtain C < 5 x 10 5 cells/mL According to this rough estimate, 
the density used in the Monte Carlo simulation (C = 10 5 cells/ml) should not lead to 
significant dipole-dipole interactions. The limit of this simple description is that it 
implicitly assumes a uniform distribution of cells at the bottom of the chamber, i.e. it 
neglects the fact that the cells, due to the DEP field, will concentrate at the electrodes 
edges (see Figure 4). In order to qualitatively take it into account this effect, it is possible 
to correct the formula using the following parameter: 



r DEP = V£ 2 



\Max 



/<V£ 2 ) 



(12) 



where VE 2 \j\4 ax and <V£ 2 > are, the highest and average value of the electric field gradient 
calculated at the bottom of the chamber. Since the dielectrophoretic force is proportional 
to the electric field gradient, r DEP approximately represents the concentration factor 
generated by the DEP force. So that we improve the previous Eq. 11 as: 



Cthreashold = 4 X 10 



hd r DEP 



cells / ml 



(13) 



i.e. the highest the concentration ratio, the lower the threshold will be. Another 
assumption of Eq. 11 is that all the cells settle at the bottom of the chamber. This is 
not generally true (see Figure 3) and it depends on the allowed deposition time and on 
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the electric fields. To take this into account we need to substitute, in Eq. 13 the chamber 
height h with the effective capturing region (z capture ): 

Cthreashoid = 4 x 10 8 ^ cells / ml 

Z capture^ ^DEP 

Zcapture Zcapture (14) 

/ ( Vz (z))dz= f {v ; Ueffiz)) d Z =t dep 

J J 6TTii m r ceU 

o o 



Where (V 1 z ll e j^z)) is the average DEP force in the vertical direction at distance z from 
the bottom of the chamber and u m is medium viscosity. In the case of the polynomial 
electrodes used and for a deposition time of 180 sec, r DEP = 85 and z cap = 280/um = (l/5)h 
so that the improved concentration threshold, to avoid dipole-dipole interaction, should 
be below 3 x 10 4 cells/ml. We performed Monte Carlo simulations in the 10 4 cells/ml 
range finding no significant evidence of cell chains formation (not shown), thus 
confirming that Eq. 14, together with Eq. 12, represent a better qualitative threshold 
to avoid cell-cell interaction requiring only a knowledge of the electric field in the 
DEP device. 



Conclusions 

In conclusion, we have demonstrated that the effects of particle-particle interactions play 
a crucial role in the kinetic evolution of colloidal systems in DEP devices even at low 
density regimes (10 5 cells/ml), lower than the ones currently used in DEP devices [30]. 

Monte Carlo methods allow for the simulation of sufficiently large systems in terms 
of size and number of particles (i.e. within the experimental scopes). The discrete 
approach (i.e. the particle resolution), as opposed to the fluid-flow methodologies, is 
the key ingredient of the method improvement. In the case of MDA-MB-231 tumor 
cells suspended in a static, low conductive, fluid under the action of a positive-DEP 
field generated by a polynomial schema we have elucidated the crucial role of particle- 
particle interactions on the trapping efficiency of the device, on the organization of cells 
in ordered chains and on the overall cell space distribution. We have also deduced a 
new qualitative concentration threshold to avoid cell-cell interaction which requires 
only a knowledge of the electric field in the DEP device. Clearly, to have an exact deter- 
mination of the concentration threshold for the specific DEP device used, MC-P kinetic 
simulations varying the cells density, such as the one proposed in this paper, need to be 
performed. 

Future works will be devoted to generalize the formalism here presented in order to 
include second-order effects such as cell-sedimentation and cell-stitching or including 
hydrodynamic forces to simulate cells distributions in dynamic separation systems. 
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